program constrained_relaxation;
uses common;
uses tools,inp,outp,ppot,decor,etools,relax;
Type
   r_xyz  = rsites3;
   r_e	  = rsites;
   r_tp	  = array[1..nsites] of integer;
Var
   i,j,k,l,m,n,o,dbg : integer;
   txt		     : text;
   r		     : r_xyz;
   ip		     : array[1..nsites] of atomrec;
   rb		     : array[1..nsites] of v33;
   itp,ito	     : r_tp;
   nsteps,omax,tinp  : integer;
   xgamma,dr2,drmax  : myreal;
   ueff,elas	     : myreal;
   uo,uo2	     : array[1..t_orb,1..atps] of myreal;
   tname,ss	     : mystring;
   rbt		     : v33;
   w1,w2,w3	     : v3;
   nofr		     : integer;
   dr_max,rcplus,cmp : myreal;
   po_max	     : i3;
   se,sel	     : rsites;
   quiet,nat,rtvl    : integer;

begin
   export_occup:=false;
   for o:=1 to t_orb do o_rel[o]:=false;
   reset(txt,'INPUT/objects.def');
   readln(txt,deconame);
   readln(txt,o); {nu_offset}
   readln(txt,dpmodel);
   readln(txt,o); if(o>0) then repre5:=true else repre5:=false;
   close(txt);
   read(nsteps,quiet,cueff,xgamma,drmax,tinp,rtvl,nofr);
   read_tvl:=(rtvl=1);
   if (nofr>0) then for o:=1 to nofr do begin
      read(i); o_rel[i]:=true;
   end;
   readln;
   readln(tname);
   {writeln('EPC_INPUT: ',nsteps:1,' ',xgamma:1:4,' ',drmax:1:4,' ',tinp:1,' ',nofr:1,' ',tname);}
   if(quiet=0) then write('RELAXATION: ');
   if (nofr=0) then begin
      if (quiet=0)then writeln(' constrained, all orbits.');
      for o:=1 to t_orb do o_rel[o]:=true;
   end else begin
      if (quiet=0) then begin
	 write(' constrained, only orbits: ');
	 for o:=1 to t_orb do if (o_rel[o]) then write(o:1,' ');
	 writeln;
      end;
   end;
   dbg:=0;
   basis;
   reset(txt,tname);
   IMPORT_TILING_BODY(tinp,txt);
   TABLE_INIT;
   if (tinp>0) then IMPORT_CONFIG(txt);
   close(txt);
   UNIT_CELL(multi);  
   for i:=1 to 3 do begin
      for o:=1 to 3 do blt[o,i]:=bl3[i,o];
   end;
   PP_INP;
   if (tinp=0) then begin
      i:=0;
      READ_DECO(i); {reads objects.def, ap.def}
      for i:=1 to tnodes do nan[i]:=0;
      for i:=1 to tnodes do NEW_NODE_OBJECTS(i,1);
   end;
   omax:=0;
   nat:=0;
   for i:=1 to tnodes do begin
      for j:=1 to nan[i] do begin nat:=nat+1;
	 for o:=1 to 3 do r[nat,o]:=bound(n_r[i,j,o]+rt[i,o]);
	 itp[nat]:=n_a[i,j,1];
	 ito[nat]:=n_a[i,j,2];
	 ip[nat]:=n_para[i,j];
	 {basis for each atom relative to object origin [node]}
	 get_atom_basis(ip[nat],rb[nat]);
	 if (ito[nat]>omax) then omax:=ito[nat];
	 if (itp[nat]=0) then nat:=nat-1 else begin end;
      end;
   end;
   if(quiet=0) then writeln('Atomic decoration: ',nat:1,' atoms.');
   OCH_STATUS;
   for i:=1 to t_orb do for j:=1 to atps do begin  {orbit energies}
      uo[i,j]:=0.0; uo2[i,j]:=0.0;
   end;
   rcplus:=1.0; {Angstrom}
   nntable(rcplus);
   cmp:=0;
   elastic_E(cmp,epot,elas,se,sel);
   writeln('E_ini=',epot/nat:1:8,' E_elas=',elas/nat:1:8,' -> estim E_min=',(epot-elas)/nat:1:8);
   if((quiet=0) and (nsteps>0)) then writeln(' step      E[eV/atom]      dR[A]');
   for i:=1 to nsteps do begin
      relax_para(1,cueff,xgamma,drmax,dr2,ueff);
      if(quiet=0) then writeln(i:5,' ',ueff/nat:15:8,' ',sqrt(dr2/nat):10:5);
   end;
   if(nsteps=0) then begin
      UPPAW(cueff,ueff,se);
   end else begin
      cmp:=0;
      elastic_E(cmp,epot,elas,se,sel);
   end;
   if (cueff>0.00001) then writeln('Ueff=',(ueff/nat):1:8,' (cueff=',cueff:1:2,')');
   if(nsteps>0) then begin
      writeln('E_rel=',epot/nat:1:8,' E_elas=',elas/nat:1:8,' -> estim E_min=',(epot-elas)/nat:1:8);
   end;
   EXPORT_XYZ_ENES('xyz.epc',se,sel);
   EXPORT_CONFIG('copa.epc');
end.
